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The infinite time-evolving block decimation (iTEBD) algorithm [Phys. Rev. Lett. 98, 070201 (2007)] allows 
to simulate unitary evolution and to compute the ground state of one-dimensional quantum lattice systems in 
the thermodynamic limit. Here we extend the algorithm to tackle a much broader class of problems, namely the 
simulation of arbitrary one-dimensional evolution operators that can be expressed as a (translationally invariant) 
tensor network. Relatedly, we also address the problem of finding the dominant eigenvalue and eigenvector of a 
one-dimensional transfer matrix that can be expressed in the same way. New applications include the simulation, 
in the thermodynamic limit, of open (i.e. master equation) dynamics and thermal states in ID quantum systems, 
as well as calculations with partition functions in 2D classical systems, on which we elaborate. The present 
extension of the algorithm also plays a prominent role in the infinite projected entangled-pair states (iPEPS) 
approach to infinite 2D quantum lattice systems. 
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I. INTRODUCTION 

The development of numerical methods to explore the prop- 
erties of strongly correlated many-body systems remains one 
of the most challenging problems in computational physics. In 
recent years, increasing attention has been paid to algorithms 
that express the state of the system as a tensor network. For 
instance, for quantum systems on a ID lattice, a matrix prod- 
uct state (MPS) [ ] represents the system's wave function in 
the density matrix renormalization group (DMRG) algorithm 
to compute ground states [ ], and in the time-evolving block 
decimation (TEBD) algorithm to simulate time evolution [ ] . 
Similarly, the tensor product state (TPS) [ ] and the projected 
entangled-pair state (PEPS) [j] have been proposed to ac- 
complish those tasks in 2D lattices, whereas the multi-scale 
entanglement renormalization ansatz (MERA) [ ] is specially 
suited to describe systems at criticality or with topological or- 
der. Finally, tensor networks can also be used to encode and 
manipulate the partition function of 2D classical lattice sys- 
tems [7, 8,9]. 

The computational cost of a simulation using tensor net- 
work algorithms is roughly proportional to the size of the lat- 
tice. However, when the system is invariant under transla- 
tions, this cost can be made independent of the system's size. 
The infinite TEBD (iTEBD) algorithm [7] exploits this fact 
to simulate unitary evolution and compute the ground state 
of a ID quantum system in the limit of an infinite lattice. 
The key idea is to encode the wave function in an infinite 
MPS (iMPS) made of a small number of tensors that are re- 
peated indefinitely and, importantly, to maintain the iMPS in 
its canonical form during the whole simulation. As a result, 
bulk properties of ID quantum systems are computed directly 
in the thermodynamic limit, circumventing more costly and 
less accurate approaches based on finite-size scaling. Other 
algorithms, such as the power wave function renormalization 
group (PWFRG) [10] or infinite DMRG (iDMRG) [11] also 
compute the ground state of infinite systems. 

A major limitation of the iTEBD algorithm is that it can 
only address unitary evolution [as explained in Sect. Ill, the 
computation of ground states with imaginary time evolution is 



a lucky exception]. Thus, the simulation of more general types 
of evolution, such as master equation evolution in a dissipa- 
tive system or imaginary time evolution to compute thermal 
states, is still restricted to finite systems [13, 14]. The reason 
lies in the fact that only unitary evolution preserves the canon- 
ical form of the iMPS. The latter is essential in order to keep 
truncation errors small during the simulation. Indeed, in the 
absence of the canonical form, truncation errors accumulate 
unnecessarily fast and ruin the simulation in places where an 
efficient iMPS description would otherwise still be feasible. 

In this paper we explain how to overcome such shortcom- 
ing. First we describe how to compute the canonical form of 
an iMPS. Then we present an extension of the iTEBD algo- 
rithm that is able to simulate a much wider class of evolution. 
Namely, it simulates the action on an iMPS of any transforma- 
tion that can be expressed as a translationally invariant tensor 
network. This includes, as particular cases, evolution in imag- 
inary time or according to a master equation. We also explain 
how to use the algorithm to compute the dominant eigenvalue 
and eigenvector [15] of any one-dimensional transfer matrix 
that decomposes as a translationally invariant tensor network. 
As an application of this, we explain how to extract correla- 
tors and local observables from the partition function of a 2D 
classical system. Finally, the extended version presented in 
this work plays a prominent role in the infinite PEPS (iPEPS) 
algorithm to simulate evolution and compute ground states in 
infinite 2D quantum lattice systems [16], as well as in certain 
implementations of MERA algorithms [6]. 

We emphasize that other algorithms can be used to address 
infinite systems, and that they are also based on or related to 
computing the dominant eigenvalue and eigenvector of a one- 
dimensional transfer matrix. This is the case, for instance, 
of the PWFRG and iDMRG algorithms [10, 1 1] to compute 
ground states in ID quantum systems and the transfer ma- 
trix renormalization group (TMRG) algorithm [ ] to evaluate 
partion functions in 2D classical systems. However, iTEBD 
differs from them at its core in two important aspects: first, 
TMRG, PWFRG and iDMRG are variational methods, while 
iTEBD amounts to a power method; second, whereas TMRG, 
PWFRG and iDMRG converge towards an infinite system by 
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adding sites to a finite lattice, in iTEBD the system is infi- 
nite from the onset. We notice that Ref. [ ], which contains 
a useful comparative study, highlights that the iTEBD is sig- 
nificantly more accurate than other proposals in determining 
ground states. Nevertheless, all these methods are of compa- 
rable interest. 

The rest of the paper is organized as follows. Sect. II ex- 
plains how to obtain the canonical form of an IMPS by or- 
thonormalizing all its bond indices. Sect. Ill presents the gen- 
eralization the iTEBD algorithm to account for non-unitary 
evolution. Sect. IV discusses an application of the algorithm 
to 2D classical lattice models and Sect. V contains some con- 
clusions. Finally, the appendix presents a detailed description 
of how to implement the algorithm for very specific forms of 
the evolution operator. 



II. CANONICAL FORM OF AN INFINITE MPS 

We consider an infinite ID lattice, where each site is labeled 
by an integer r G Z and described by a Hilbert space of fi- 
nite dimension d. The lattice is in a pure state | ^) G (8)rGZ 
that is invariant under translations by n sites (and multiples 
thereof). Following Ref. [ ], we represent |^) using an 
IMPS, which in the simplest case (n = 1) consists of a pair 
of tensors {F, A}, see Figs, (l.i)-(l.n). Here F is made of 
complex coefficients FJ^^, with two bond indices a and f3 
{a^P = 1, • • • , x) and one physical index i (i = 1, • • • , d) 
that labels an orthonormal basis in C^, whereas A is a diagonal 
matrix with non-negative diagonal elements • The integer 
X is known as the rank of the IMPS. 

We say that an IMPS {F, A} is in its canonical form [ , 
17] when, on each bond, index a is related to the Schmidt 
decomposition of |^), 



(1) 



that is, when the diagonal matrix A contains the decreasingly 
ordered Schmidt coefficients (Ai > A2 > • • • > A^ > 0) and 
a labels the Schmidt vectors, which form orthonormal sets, 
(<l>^|<l>^,) = (^^1^^,) = 6aa'' In terms of the matrices 
and L, defined as 
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Figure 1 : (color online) (i) Diagrammatic representation of the ten- 
sors F and A that form an iMPS for a pure state |^) of a transla- 
tionally invariant chain, (ii) The iMPS is actually an infinite one- 
dimensional tensor network consisting of alternating copies of the 
tensors F and A. In the diagram, a leg shared by two tensors cor- 
responds to an index over which there is an implicit sum. (Hi) Dia- 
grammatic representation of the conditions in Eqs. (4)-(5), which are 
fulfilled if and only if the iMPS {F, A} is in its canonical form, (iv) 
Tensor network corresponding to the expectation value (^|0^^^|^} 
of an operator 0^^\ acting on site r, when |^) is represented by an 
iMPS {F, A} in its canonical form. As a result of the orthonormality 
of the Schmidt vectors in Eq. (1), only a very small number of ten- 
sors need to be considered, (v) Tensor network for the computation 
of the two-point correlator (^|0^^^0^*^ |^) when |^) is represented 
by an iMPS {F, A} in the canonical form. 



where 7/ G C. In other words, the identity operator I a a' = 
(^Q,^/ is a right (left) eigenvector of matrix R (respectively L) 
with eigenvalue r], see Fig.(l.zn)- Incidentally, r] is the domi- 
nant eigenvalue [15] of both R and L, and is equal to 1 if and 
only if 1^) is normalized [l ^]. 

Two good reasons to express an iMPS in its canonical form 
are the following. First, it facilitates the computation of ex- 
pectation values for local operators. For instance, for 0^^^ an 
operator acting on site r, orthogonality of the Schmidt bases 
implies that (^|0M |^) is simply 



(^|0H|^) 



a,i,(3 



(6) 



which can be computed in 0{d^y^) time. Similarly, the ex- 
pression for a two-point correlator (^|0M0[^]|^) involves 
only of the order of I s — r I tensors, see Figs. {\.iv)-{\.v). Sec- 
ond, as we will discuss in Sect. Ill, the canonical form simpli- 
fies the truncation of bond indices, a process that is necessary 
in order to prevent the rank of the iMPS from growing during 
a simulation. 

Theorem 1 of Ref. [ ] explains how to bring an MPS to 
its canonical form in the case of a finite chain. This is done 
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by orthonormalizing bond indices, starting from one boundary 
of the chain and progressing through the whole system, with 
a cost proportional to its length. In the infinite case, we use 
translational invariance to reduce this cost to a constant. More 
specifically, given an IMPS {F, A} for |^), we can obtain a 
canonical form {F^, A'} for the same state through three steps, 
illustrated in Fig. (2): 

{i) Find the matrix Vr that is the dominant right eigenvector 
[15] of R, in the sense of Fig. (2.z), with dominant eigenvalue 
T] e C (here r] is assumed to be unique [ ]). Similarly, find 
the matrix Vl that is the dominant left eigenvector of L, which 
also has eigenvalue 77 [we use a large-scale, non-Hermitian 
eigenvalue solver [20], such as an Arnoldi method, and exploit 
the tensor network structure of R and L] . Decompose matri- 
ces Vr and Vl, which are Hermitian and non-negative (since 
they originate in the scalar product of a set of non-orthogonal 
vectors), as squares Vr = XX^ and Vl = Y^Y. For in- 
stance, if Vl = WDW^ is the eigenvalue decomposition of 
Vl, then Ft = wVD and Y = VDW"^ [ ]. 

{ii) Introduce the two resolutions of the identity matrix 
I = (r^)-iy^ and I = XX-^ in the bond indices of 
the IMPS as indicated in Fig.(2.n). Then, compute the sin- 
gular value decomposition of the product Y^XX, namely 
Y^XX = UX'V, where U^V are unitary and the diagonal 
matrix A' contains the Schmidt coefficients of |^). 

(ill) Arrange the remaining tensors V, X~^, F, (Y^)~^ 
and U into a new tensor F' as in Fig. (l.iii). 
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Figure 2: (color on-line) (i-iii) The three steps involved in the com- 
putation of the canonical form {F^ A^} for an iMPS {F, A} for |^), 
as explained in the text. In particular, A^ is obtained in step (ii) as 
the singular values of Y^ XX, whereas F^ is defined in step (Hi), 
(iv) Both {F, A} and {F^, A^} represent the same state |^), as can 
be verified by direct substitution. 



All the previous manipulations can be implemented with 
computational cost scaling as 0{dx^)- A proof that the re- 
sulting IMPS {F^, A'} is indeed in the canonical form is given 
in Fig. (3). 

We can now analyze the case where |^) is invariant under 
translations by n > 1 sites. For n = 2, the state |^) is rep- 
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Figure 3: (color online) Proof that {F^A^}, the iMPS obtained from 
{F, A} by re-orthonormalizing its bond indices following Fig. (2), is 
indeed in the canonical form. The diagram only shows the condition 
of Eq. (4) for B! . The proof of Eq. (5) for L' is analogous. 



resented by an iMPS that consists of four alternating tensors 
{F^, A^, F^, A^}, where A and B denote odd and even sites 
in the chain [12]. The canonical form {F^', A^', F^', A^'}, 
defined as before to correspond to the Schmidt decomposition 
at each bond, can be obtained as follows. First we coarse- 
grain the chain by regarding each pair of sites AB as a single 
site and represent |^) with an iMPS (F, A) as in the n = 1 
case. Then we transform the coarse-grained iMPS (F, A) into 
its canonical form (F^, A'). Finally, we split F' into three ten- 
sors r-^ , and by means of a singular value decom- 
position, see Fig (4). These steps can be implemented with 
a computational cost that scales as 0((i^x^)- The case of a 
generic n is addressed similarly, and the computational cost 
scales as 0{n^x'). 
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Figure 4: (color online) Construction of the canonical form {F"^ , 
A^', F^', A^'} from an iMPS {F^, A^, F^, A^} that is invariant 
under translations by n = 2 sites: (i) Coarse-grained iMPS {F, A}. 
{ii) Canonical form {F^,A^} for the coarse-grained iMPS{F,A}. 
(iii) X^ corresponds to A^ (iv) X^ is obtained through a singular 
value decomposition, from where also {v) F^ and F^ are obtained 
after minor manipulations. 
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III. SIMULATION OF NON-UNITARY EVOLUTION 

In this section we discuss how to update the iMPS for state 
I ^) after a gate G acts on the entire lattice. That is, we aim to 
build an IMPS for the resulting state | = G| ^) . We assume 
that G is expressed as a one-dimensional tensor network (of 
some sort) that is invariant under translations by n sites, see 
Fig. (5) for several examples. As a remark, let us mention 
that non-unitary gates such as the ones in Fig. (5) appear in 
the so-called iPEPS algorithm to simulate 2D quantum lattice 
systems [ib]. 
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Figure 5: (color online) Four types of gates G acting on an iMPS: {%) 
infinite matrix product operator (iMPO); (ii) product of two-site op- 
erators; [iii) product of two-to-one- site operators; {iv) Tensor prod- 
uct of one-to-two- sites operators. Detailed explanations on how to 
update the iMPS for cases (ii)-(iv) can be found in the appendix. 
Notice that gates (Hi) and (iv) change the number of sites in the lat- 
tice. These gates appear e.g. when coarse-graining (or fine-graining) 
the chain. 

We focus again on the case n = 1 and, for concreteness, we 
assume G is specified by an infinite matrix product operator 
(iMPO) as in Fig. (5.i). This iMPO is represented by a tensor 
a of complex components a^^, where i and j are physical 
indices and fi and v (fi^v = 1, • • • , /^) are bond indices. The 
update occurs in three steps, illustrated in Fig. (6): 

(I) Contraction: the tensors {F, A} for |^) are contracted 
with the tensors that specify the gate G, producing an iMPS 
{f,A}for|*')> 



a/3 



(7) 



Here indices a and (a, = 1, • • • , x) are defined as a = 
(a, /i) and (3 = (/5, u). Notice that the rank x of the new 
iMPS, X = ^X' is larger than the rank x of the initial iMPS. 
The computational cost of this step is 0{dP' i<i^x^)- 

(II) Orthogonalization: the iMPS {f, A} for |^') is brought 
into its canonical form {F', A'} with a cost that scales as 

(III) Truncation: a final iMPS {F', A'} is obtained from 
the canonical form {F', A'} by truncating all bond indices. In 



particular, on each bond we preserve the first % values of the 
index, corresponding to the x largest Schmidt coefficients. 

ATA 

x^^ x^^x x'^^x x'^x 



Figure 6: (color online) Sequence of transformations (I) -(III) that 
produce a truncated iMPS {F'A'} for |^') = Gl^) from an iMPS 
{r, A} for 1^), as explained in the text. Notice in particular that, as 
a result of the truncation step, the initial and final iMPSs have the 
same rank x- 

The net result is an approximate iMPS {F^, A'} for |^'), 
obtained with a total computational cost 0{d?H?x^ + dn^x^) 
[22]. The truncation step is necessary in order to keep the 
rank x (and therefore the computational cost) constant during 
a simulation, where typically not just one gate G but rather 
of a whole series {Gi, G2, • ' ' } is sequentially applied to the 
chain. The truncation of the bond indices introduces an error 
that is hard to evaluate in an infinite system. Here we apply, 
to all bond indices simultaneously, the truncation scheme that 
is known to be optimal when applied only to one bond index 
[23]. 

For n > 1, as is the case e.g. in Fig. (5.u), we can again 
coarse-grain the system and proceed as in the n = 1 case, 
which will result in a iMPS {F', A'}. Then F' is broken into 
several other tensors {T'^\ X"^\t^\ ■ ■ process that may 
require additional truncations. See the appendix for a detailed 
analysis of some particular cases. 

We are therefore able to address non-unitary evolution on 
an infinite chain. When the gate G breaks into a row of 
two-site gates as in Fig. (5.ii), and each two-site gate is uni- 
tary, then the canonical form of the iMPS is preserved (up to 
truncation errors) without need of the orthogonalization step, 
recovering the original formulation of the iTEBD algorithm 
[ ]. Notice that in Ref. [ ] the algorithm is also used to 
compute the ground state of the system. This is done by sim- 
ulating (non-unitary) imaginary time evolution 



l*r) = 



exp(-ifr)l^o) 
|exp(-FT)|4'o) 



(8) 



where H is the Hamiltonian of the infinite chain and |^o) 
some initial state, and by exploiting the fact that under proper 
circumstances the ground state |^) of is the fixed point of 
such evolution. 



1^) = lim 



exp(-ij^r)l^o) 
|exp(-i7r)|^o)| 



(9) 



We emphasize that such calculation succeeds thanks to a 
fortunate combination of favorable, unlikely circumstances. 
When the simulation is performed using small time steps, two- 
site gates that are close to the identity operator are used. These 
gates destroy the canonical form of an initial iMPS, but they 
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leave its bond indices in a (non-orthonormal) basis that still 
seems to lead to reasonably small errors during their trunca- 
tion. One would expect the bond bases to become less and 
less adequate for truncation over time, as the accumulated r 
increases, since the overall evolution exp(— i^r) departs more 
and more from the identity. But it turns out that the singu- 
lar value decomposition used in order to update the IMPS at 
each time step has the effect of reorganizing the indices fa- 
vorably, partially compensating the non-unitary effects [24]. 
Finally, all the excessive truncation errors introduced during 
the simulation are washed away at its final stages, where in- 
creasingly small time steps are used. These have the intended 
effect of reducing Suzuki-Trotter errors [ ], but they also im- 
ply that the gates become almost unitary (that is, very close to 
the identity). One can see that, as a result, by the end of the 
simulation the IMPS approximation for the ground state is not 
only accurate, but it is also very close to the canonical form. 

Thanks to actively transforming the IMPS into its canon- 
ical form, the present extension is not restricted to unitary 
evolution and can be applied to a wider range of ID prob- 
lems. In particular, it can be used to simulate master equation 
evolution and to compute thermal states using the mixed state 
formalism of Refs. [13]. Importantly, it can also be used to 
manipulate the state of an infinite 2D lattice, both for classi- 
cal (see example below) and quantum systems [ ]. This is 
achieved after the 2D problem is recast into that of finding the 
dominant eigenvalue 9 and dominant eigenvector [ ] |^) of 
a ID transfer matrix T that decomposes into a finite sequence 
{Gi, G2, • • • , Gm} of gates. The dominant eigenvector satis- 
fies 



The system's partition function reads 



1^) = lim 

p^oo 



|T^|^o)| 



(10) 



and is obtained by simulating the repeated application of T on 
an initial state |^o), until converge is attained. The dominant 
eigenvalue can be obtained from the dominant eigenvector | ^) 
and any other vector |<l>), (^1^)7^0, since 



(<I>|T|^) 



(11) 



In the next section we provide an explicit example of cal- 
culation of dominant eigenvalue and eigenvector of a one- 
dimensional transfer matrix. 



IV. EXAMPLE: 2D CLASSICAL SYSTEMS 

In this section we explain how the above algorithm can be 
used to compute the partition function, local observables and 
two-point correlators of a classical spin system. We consider 
an infinite 2D lattice where each site, labeled by a vector r, 
contains a d-dimensional spin that interacts with nearest 
neighbor spins according to a Hamiltonian K, 



K{{s}) 



{r^r') 



1). 



(12) 



Z(/3) = ^e- 



■pK{{s}) 



En 

{s} {r,r') 



g-/3X.(.l^al^-'l)^ (13) 



where f3 is the inverse temperature. For concreteness, we con- 
sider a square lattice with an isotropic interaction K2. Let 
^/Q denote the squared root of the Hermitian matrix Qss' = 
ex.p{—PK2{s^ s')) [25]. We can express the partition function 
as the contraction of an infinite 2D tensor network specified 
by a single tensor a, 

aijkl = ^{VQ)is{VQ)js{VQ)qs{VQ)ls (14) 



that is repeated on all sites [ ], see Fig. (7). The above tensor 
can be computed in 0{d^) time. 



3^ 



Figure 7: (color online) The partition function of a 2D classical 
system can be written as the contraction of a 2D tensor network. In 
the case of an infinite square lattice with isotropic and homogeneous 
interactions, this tensor network consists of infinitely many copies of 
the tensor a in Eq. (14). 

We now introduce an infinite ID transfer matrix T consist- 
ing of one row of tensors a, see Fig. (8). Then we have 



Z{p) = lim tr(T^) = lim / 

p^oo p^oo 



(15) 



where 9 is the dominant eigenvalue of T. Let \^u) and \^d) 
be the corresponding (up and down) eigenvectors, 

T\<Su) = e\<Su), {^D\T = e{^Dl (16) 

that we normalize to {^d\^u) = 1- Then 

e = {^d\T\^u) = triW) = lim w", (17) 

q^oo 

where co is the dominant eigenvalue of matrix W defined in 
Fig. (8), and we finally have 



lim 

p,q^oo 



(18) 



Therefore, in order to evaluate the partition function Z{f3), 
we will first construct an IMPS {F^, A^} for and an 

IMPS {r^, A^} for I^^d) by iteratively applying the transfer 
matrix T on an initial state |^o) (cf. Eq. (10)). Specifically, 
we use the iTEBD algorithm as discussed in the previous sec- 
tion to simulate the state 



0/ 



||Tf|*o)||' 



p= 1,2,-- 



(19) 
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for increasing values of p, until the resulting IMPS has con- 
verged within some agreed precision [ ]. The computation 
is approximate, in that an IMPS with finite rank x will be used 
to represent the dominant eigenvectors, which in general may 
only be represented exactly with an infinite rank x- Notice 
that if K2 is isotropic, the transfer matrix can be made Her- 
mitian, in which case = \^u) [otherwise I^d) also 

needs to be computed]. From the converged iMPSs {F^, A^} 
and {F^, A^} we can construct matrix W. The dominant 
eigenvalue uoofW can then be computed using a large-scale 
eigenvalue solver and exploiting its tensor network structure 
in 0{d'^X^ ^ d'^X^) time. 



m m 




Figure 8: (color online) Computation of the partition function 
through Eq. (18). (i) is the up dominant eigenvector of the 

transfer matrix T, with dominant eigenvalue 0. (ii) Similarly, |^d) 
is a down dominant eigenvector of T, with dominant eigenvalue 6. 
(Hi) \(Pr) is the right dominant eigenvector of matrix with dom- 
inant eigenvalue to. (iv) \(Pl) is the left dominant eigenvector of 
matrix W, with dominant eigenvalue u. 

On the other hand, for any function f{s) of one spin, the 
expectation value 

is, up to the factor 1/Z{(3), also given by the contraction of 
an infinite 2D tensor network, obtained from that for by 
replacing tensor a on site r with tensor b, 

hjki = ^/(5)(\/Q)^.(yQ),.(v/Q)g.(yQ)^. , (21) 




Figure 9: (color online) The expectation value {f{s^ )) from Eq. 
(20) is the ratio of the contraction of two infinite 2D tensor networks. 
By introducing the dominant eigenvectors of the one-dimensional 
transfer matrix T, we can rewrite (/(s^^)) as the ratio of the trace 
of two infinite ID tensor networks. Finally, by introducing the dom- 
inant eigenvectors of matrix W , we obtain a ration of two simple 
tensor networks. 

and of the matrix W. We have already indicated how to 
proceed in the computation of these quantities. 

Similarly, we can build a tensor network for the expectation 
value of the correlator 

{.} 

(22) 

by replacing the tensor a in sites r and r' with appropriate ten- 
sors h and h' and proceeding in a similar way as the previous 
case, see Fig. (10). Notice that we assume that sites r and r' 
lie on the same row of the lattice. 




again computable in 0{(fi) time. As illustrated in Fig. (9), 
is eventually written as the ratio of two small ten- 
sor networks. These tensor networks are expressed entirely 
in terms of: tensors a and h\ tensors {F^, A^} and {F^A^} 
defining the dominant eigenvectors \ ^u) and | ^i)) of the one- 
dimensional transfer matrix T; and the dominant vectors | Lpji) 



Figure 10: (color online) Following steps analogous as those of 
Fig. (9) for the expectation value (/(s^^)), the two-point correla- 
tor (/(s^^)5'(s^^ 0) can also be reduced to the ratio of two simple 
tensor networks. 

Fig. (11) shows the magnetization per site m = {s^'^^) for 
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Figure 11: (color online) Magnetization per lattice site for the 
infinite-size 2D classical Ising model at different temperatures f3. 
The exact solution m = (1 — ((sinh (2/3))"^))^^^ has been included. 
The numerical results have been obtained by approximating the dom- 
inant eigenvectors of the one-dimensional transfer matrix T with an 
IMPS of rank = 40. The inset shows the relative error. 
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Figure 12: (color online) Two-point correlators for the infinite- size 
2D classical Ising model at critical temperature, Pc — ^ log(l+A/2), 
along a row of the square lattice. The exact solution, which scales as 
(^g[r\ g[r']^^ ^ c|r — r'l"^/^, has been included. The numerical 
results have been obtained by approximating the dominant eigenvec- 
tors of the one-dimensional transfer matrix T with an IMPS of rank 
= 40, X = 60 and x = 80. 



V. CONCLUSIONS 

In this paper we have explained how to extend the iTEBD 
algorithm so that it can be applied to simulate any evolution 
that can be expressed as a sequence of one-dimensional tensor 
networks. The key new ingredient is a recipe to rewrite any 
IMPS in the canonical form, as required in order to properly 
truncate the bond indices. 

The iTEBD algorithm can therefore be applied to simulate 
not only unitary evolution, but also master equation evolu- 
tion and imaginary time evolution. It can also be used to find 
the dominant eigenvalue and dominant eigenvector of a one- 
dimensional transfer matrix. This last application is particu- 
larly relevant in order to analyze the partition function of 2D 
classical models, as we explained, and it also plays a promi- 
nent role in the iPEPS algorithm for 2D quantum systems [16] 
and some of the MERA algorithms [ ]. 
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Appendix A: ALGORITHMS FOR THE EVOLUTIONS IN 
FIG.(5) 

In this appendix we explain in some detail how to im- 
plement the iTEBD algorithm for some particularly relevant 
choices of the gate G. In Sect. Ill we analyzed the case where 
the system is invariant under translations by n = 1 lattice 
sites. Here we consider the case n = 2 for the four different 
types of gate represented in Fig. (5). Most of the information 
contained in the appendix can already be derived from the re- 
sults of Sects. II and III, but we write it explicitly for each 
gate for the sake of clarity. 



1. Matrix product operator, Fig.(5.i) 



the 2D Ising model, defined by 



K2{s,s') = -s 5,s' = ±l, (23) 



at different values of p. We have used an IMPS of rank x = 40 
to represent the dominant eigenvectors and and 
proceeded as explained above. It is noteworthy that the nu- 
merical results reproduce the exact behaviour of m with small 
relative error. Furthermore, Fig. (12) shows the decay with 
\f — f*\ of the spin-spin correlator (s^^st^]) at the critical 
point, (3c = ^ log(l + V^)- In this case we have used an 
iMPS of rank x — 40, 60, 80. Remarkably, the numerical re- 
sults reproduce the correct power-law decay ~ | r — | ~ for 
distances up to thousands of spins, with increasing accuracy 
as X increases. 



Let us assume that the iMPS for |^) is given by some ten- 
sors r^, for odd sites and tensors F^, for even sites. 
We consider the evolution under a MPO invariant under trans- 
lations of two chain sites, 

G= E --Mr'^r'C''--, (Ai) 

Q!,/3,7,5,p,... 

where for each value of a, /3 = 1, 2, . . . , the one- site oper- 
ators a[^j3 : and : H^^+i] ^ ^[^+^1 act on 
sites i and z + 1, see Fig.(5.z). The updating procedure of the 
iMPS can be expressed as follows: 

{i) Compute tensor 6i as indicated in Fig.(13.z), with bond 
dimension hiX- 

{a) Find the matrix Vr that is the right dominant [15] 
eigenvector of R (in the sense of Fig. (13. n)) with dominant 
eigenvalue r] e C (assumed to be unique) [19], where R is 
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obtained by contracting 61 with its complex conjugate 9* 
as shown in Fig. (13. n) [use large-scale, non-Hermitian eigen- 
value solver, such as Arnoldi methods, and exploit the tensor 
network structure of R] . Then, decompose matrix Vr (which 
is Hermitian and non-negative) as the square Vr = XX^. For 
instance, if Vr = WDW^ is the eigenvalue decomposition of 
Vr, then X = W^. 

{Hi) Compute tensor 62 as indicated in Fig.(13.zn). with 
bond dimension kx- 

{iv) Find the matrix Vl that is the left dominant [ l^] eigen- 
vector of L (in the sense of Fig. (13. zv)) with dominant eigen- 
value r G C (assumed to be unique) [19], where L is obtained 
by contracting 62 with its complex conjugate 6 J as shown in 
Fig.(13.zv) [use large-scale, non-Hermitian eigenvalue solver, 
such as Arnoldi methods, and exploit the tensor network struc- 
ture of L]. Then, decompose matrix Vl (which is Hermitian 
and non-negative) as the square Vl = Y^Y, in the same way 
as was done for matrix Vr. 

{v) Compute tensor O as indicated in Fig.(13.v), with bond 
dimension /^x- 

(vi) Introduce the two resolutions of the identity matrix 
I = {Y^)-^Y^ and I = XX'^ in the bond indices of 
tensor 9 as indicated in Fig.(13.vi). Then, compute the sin- 
gular value decomposition Y^ X = UX^'V, leading to new 
Schmidt coefficients Truncate these new Schmidt coeffi- 
cients by keeping only the x largest ones, and normalize them 
so that the sum of their squared values is 1 . 

{vii) Compute tensor E as indicated in Fig.(13.vn). 

{iix) Group the indices of S according to a single index for 
the left-hand side and a single index for the right-hand side, 
and compute the singular value decomposition as indicated in 
Fig.(13.nx). This leads to two isometric tensors P and Q, and 
new Schmidt coefficients Truncate these new Schmidt 
coefficients by keeping only the x largest ones, and normalize 
them so that the sum of their squared values is 1 . 

{ix) Obtain new matrices F^' and F^' as indicated in 
Fig.(13.ix). 

The above sequence of steps has a computational cost of 
0{(fhi^X^) ill t™^- 



2. Tensor product of two-site operators, Fig.(5.ii) 

Consider an IMPS for state | ^) with bond dimension x that 
is invariant under shifts of two chain sites. This IMPS is then 
defined by tensors F^, for odd sites and tensors F^ , A^ for 
even sites. The evolution operator that we consider is given by 
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where a[^'^+i] : H^^ H^'"'^^ is a two-body 

operator acting on the two contiguous sites i and i + 1 of the 
IMPS, see Fig.(5.n). The algorithm to update the IMPS is as 
follows: 

(i) Compute tensor 61 as indicated in Fig.(14.z), with bond 
dimension x- 
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Figure 13: (color online) Steps in the updating of the iMPS after the 
action of a MPO. 



(ii) Find the matrix Vr that is the right dominant [15] 
eigenvector of R (in the sense of Fig.(14.n)) with dominant 
eigenvalue r] e C (assumed to be unique) [ ], where R is 
obtained by contracting 9i with its complex conjugate 9^ 
as shown in Fig.(14.n) [use large-scale, non-Hermitian eigen- 
value solver, such as Arnoldi methods, and exploit the tensor 
network structure of R] . Then, decompose matrix Vr (which 
is Hermitian and non-negative) as the square Vr = XX^. For 
instance, if Vr = WDW^ is the eigenvalue decomposition of 
VR,±enX = W^/D. 

{Hi) Compute tensor 62 as indicated in ¥ig.(l4.iii), with 
bond dimension x- 

{iv) Find the matrix Vl that is the left dominant [15] eigen- 
vector of L (in the sense of Fig.(14.iv)) with dominant eigen- 
value r G C (assumed to be unique) [ ], where L is obtained 
by contracting 62 with its complex conjugate 6 J as shown in 
Fig. (14.2^) [use large-scale, non-Hermitian eigenvalue solver, 
such as Arnoldi methods, and exploit the tensor network struc- 
ture of L] . Then, decompose matrix Vl (which is Hermitian 
and non-negative) as the square Vl = Y^Y, in the same way 
as was done for matrix Vr. 
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(v) Compute tensor 6 as indicated in Fig.(14.v), with bond 
dimension x- 

(vi) Introduce the two resolutions of the identity matrix 
I = {Y^)-^Y^ and I = XX'^ in the bond indices of 
tensor 6 as indicated in ¥\g.{\A.vi). Then, compute the sin- 
gular value decomposition X = UX^'V, leading to new 
Schmidt coefficients Truncate these new Schmidt coeffi- 
cients by keeping only the x largest ones, and normalize them 
so that the sum of their squared values is 1 . 

{vii) Compute tensor E as indicated in ¥ig.(\A.vii). 

{iix) Group the indices of E according to a single index for 
the left-hand side and a single index for the right-hand side, 
and compute the singular value decomposition as indicated in 
¥ig.{\A.iix). This leads to two isometric tensors P and Q, and 
new Schmidt coefficients Truncate these new Schmidt 
coefficients by keeping only the x largest ones, and normalize 
them so that the sum of their squared values is 1 . 

{ix) Obtain new matrices T^' and T^' as indicated in 
Fig.(14.ix). 

The computational cost of the above sequence of steps is 
0{dt^X^). Also, and as expected, if a[*'*+^] is a unitary oper- 
ator then this procedure corresponds exactly to the updating 
rules of the standard iTEBD algorithm. 
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3. Tensor product of two-to-one-site operators, Fig.(5in) 

Our concern now is the evolution of an IMPS under an op- 
erator G that is the tensor product of two-to-one- site operators 



(iix) 



G 



(A3) 



ie odd 



where a[^'^+^] : H^^ H^'^^^ n^i^^^^/'^] is a two-to-one- 
site operator acting on two contiguous sites i and i + 1 of the 
IMPS, and which maps the two sites to a new site (z + l)/2, 
see Fig.(5in)- Again we assume that the IMPS is defined by 
r^, for odd sites and , A^ for even sites. The algorithm 
to update the IMPS is as follows: 

{i) Compute tensor 6i as indicated in Fig.(15.i), with bond 
dimension x- 

(a) Find the matrix Vr that is the right dominant [1^] 
eigenvector of R (in the sense of Fig.(15.n)) with dominant 
eigenvalue r] e C (assumed to be unique) [19], where R is 
obtained by contracting Bi with its complex conjugate 6^ 
as shown in Fig. (15. n) [use large-scale, non-Hermitian eigen- 
value solver, such as Arnoldi methods, and exploit the tensor 
network structure of R] . Then, decompose matrix Vr (which 
is Hermitian and non-negative) as the square Vr = XX^. For 
instance, if Vr = WDW^ is the eigenvalue decomposition of 
Vr, then X = WVD. 

{Hi) Compute tensor 02 as indicated in Fig. (15. in), with 
bond dimension x- 

(iv) Find the matrix Vl that is the left dominant [ ] eigen- 
vector of L (in the sense of Fig. (15. iv)) with dominant eigen- 
value r G C (assumed to be unique) [ ], where L is obtained 
by contracting 62 with its complex conjugate 9 J as shown in 
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Figure 14: (color online) Steps in the updating of the iMPS after the 
action of the tensor product of two- site operators. 



Fig.(15.zv) [use large-scale, non-Hermitian eigenvalue solver, 
such as Arnoldi methods, and exploit the tensor network struc- 
ture of L] . Then, decompose matrix Vl (which is Hermitian 
and non-negative) as the square Vl = Y'^Y, in the same way 
as was done for matrix Vr. 

{v) Compute tensor 6 as indicated in Fig.(15.v), with bond 
dimension x- 

{vi) Introduce the two resolutions of the identity matrix 
I = {Y^)-^Y^ and I = XX in the bond indices of ten- 
sor 6 as indicated in ¥ig.(l5.vi). Then, compute the singular 
value decomposition Y^X = UyV, leading to new Schmidt 
coefficients A^ Truncate these new Schmidt coefficients by 
keeping only the x largest ones, and normalize them so that 
the sum of their squared values is 1 . 

(ix) Obtain a new matrix as indicated in Fig.(15.vn). 

The above procedure has a computational cost of 0{d^x^)- 
In the end, the action of the two-to-one- site gates a^^^ on the 
IMPS can be computed exactly without further truncation of 
the bond indices, and is such that the obtained IMPS for the 
evolved state 1^0 is invariant under translations of one chain 
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Figure 15: (color online) Steps in the updating of the iMPS after the 
action of the tensor product of two-to-one- site gates. 



site instead of two. 
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4. Tensor product of one-to-two-sites operators, Fig.(5.iv) 

Contrary to the previous cases, we consider now the situa- 
tion in which the iMPS for state | ^) is defined by one tensor T 
and one Schmidt vector A, so that it is invariant under shifts of 
one chain site. At this point we wish to update the iMPS after 
the action of a tensor product of one-to-two- sites operators G 

G = (g)a[^] , (A4) 

i 

where a^^ : H^^ H^^^"^] (g) H^'^''^ is a one-to-two-sites 
operator acting on one site i of the iMPS, and which maps the 
site to two new sites 2i — l and 2i, see ¥ig.(5.iv). The steps to 
follow to update the iMPS are: 

{i) Compute tensor 6i as indicated in Fig.(16.z), with bond 
dimension x- 

(ii) Find the matrix Vr that is the right dominant [15] 
eigenvector of R (in the sense of Fig.(16.n)) with dominant 
eigenvalue r] e C (assumed to be unique) [ ], where R is 
obtained by contracting 6i with its complex conjugate 
as shown in Fig.(16.n) [use large-scale, non-Hermitian eigen- 
value solver, such as Arnoldi methods, and exploit the tensor 
network structure of R] . Then, decompose matrix Vr (which 
is Hermitian and non-negative) as the square Vr = XX^. For 



Figure 16: (color online) Steps in the updating of the iMPS after the 
action of the tensor product of one-to-two- sites gates. 



instance, if Vr = WDW^ is the eigenvalue decomposition of 
Vr, then X = WVD. 

{Hi) Compute tensor 62 as indicated in Fig.(16.zn). with 
bond dimension x- 

{iv) Find the matrix Vl that is the left dominant [15] eigen- 
vector of L (in the sense of Fig.(16.zv)) with dominant eigen- 
value r G C (assumed to be unique) [ ], where L is obtained 
by contracting 62 with its complex conjugate 9 J as shown in 
¥\g.{\6.iv) [use large-scale, non-Hermitian eigenvalue solver, 
such as Arnoldi methods, and exploit the tensor network struc- 
ture of L] . Then, decompose matrix Vl (which is Hermitian 
and non-negative) as the square Vl = Y^Y, in the same way 
as was done for matrix Vr. 

{v) Compute tensor 6 as indicated in Fig.(16.v), with bond 
dimension x- 

{vi) Introduce the two resolutions of the identity matrix 
I = {Y^)-^Y^ and I = XX'^ in the bond indices of 
tensor 9 as indicated in ¥ig.(l6.vi). Then, compute the sin- 
gular value decomposition Y^X = UX^'V, leading to new 
Schmidt coefficients Truncate these new Schmidt coeffi- 
cients by keeping only the x largest ones, and normalize them 
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so that the sum of their squared values is 1 . 

{vii) Compute tensor S as indicated in Fig.(16.vn). 

(iix) Group the indices of S according to a single index for 
the left-hand side and a single index for the right-hand side, 
and compute the singular value decomposition as indicated in 
¥ig.(l6.iix). This leads to isometric two tensors P and Q, and 
new Schmidt coefficients Truncate these new Schmidt 
coefficients by keeping only the x largest ones, and normalize 
them so that the sum of their squared values is 1 . 



(ix) Obtain new matrices T^' and T^' as indicated in 
Fig.(16ix). 

The computational cost of the above steps is 0{d^x^)- Sim- 
ilarly to the case of the previous section, the translational in- 
variance of the original IMPS for |^) has been modified, in 
a way that the obtained IMPS representation for the evolved 
state 1^') has periodicity under shifts of two chain sites in- 
stead of one. 
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